Autowaves in the model of avascular tumour growth 



O 
O 
(N 



oo 

h; 

d ; 
cr: 



> 

o 
o 



- 

X 



A.V. KolobovQ V.V. GubernovQ and A.A. Polezhae\|l 
P.N. Lebedev Physical Institute of Russian Academy of Sciences, 53 Leninsky prospect Moscow, 119991 RUSSIA 

(Dated: February 1, 2008) 

A mathematical model of infiltrative tumour growth taking into account cell proliferation, death 
and motility is considered. The model is formulated in terms of local cell density and nutrient (oxy- 
gen) concentration. In the model the rate of cell death depends on the local nutrient level. Thus 
heterogeneous nutrient distribution in tissue affects tumour structure and development. The exis- 
tence of automodel solutions is demonstrated and their properties are investigated. The results are 
compared to the properties of the Kolmogorov-Petrovskii-Piskunov and Fisher equations. Influence 
of the nutrient distribution on the autowave speed selection as well as on the relaxation to automodel 
solution is demonstrated. The model adequately describes the data, observed in experiments. 



I. INTRODUCTION 

According to experimental data tumour growth can 
be naturally subdivided into two stages The initial 
stage, called avascular growth, is characterized by the tu- 
mour consumption of crucial nutrients, such as oxygen or 
glucose, basically via diffusion. Depending on nutrients 
concentration tumour cells are supposed to be in one of 
the three states: proliferating, resting or dead. An avas- 
cular tumour is a compact spherically symmetric colony 
of cancer cells with the necrotic region in the center. Such 
a rigid structure holds until the tumour size does not ex- 
ceed several millimeters. Later in response to chemical 
signals from cancer cells new capillary blood vessels start 
to grow to provide better nutrients supply. This process, 
called angiogenesis, indicates the second phase of tumour 
development - vascular growth. 

Although a significant progress in experimental meth- 
ods has been achieved by now, they still cannot give a 
conceptual framework within which all the existing data 
of the tumour development can be fitted 0. Thereupon 
the interest in mathematical modelling of the tumour 
growth and progression has rapidly grown in the last 
decades. The majority of these models consider avas- 
cular tumour growth or growth of multicellular tumour 
spheroids (MTS) which are prevalent experimental in 
vitro models of avascular tumours. 

There are several basic classes of mathematical models 
for the avascular tumour or MTS growth. One of them 
corresponds to "reaction-diffusion type" models which 
generally consist of an ordinary differential equation for 
the tumour volume, coupled to one or more parabolic 
partial-differential equations describing the distribution 
of nutrients and growth inhibitory factors within the tu- 
mour d, 0, [1| . In models of this type nutrient concen- 
tration provides a mitotic and/or death index and the 
equation for the tumour volume follows from the mass 
conservation law. Thus, reaction-diffusion models do not 
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consider any cell motion. Also these models can deal 
only with spherically or cylindrically symmetric spatial 
structure of a tumour. 

In convection-domination models a tumour is consid- 
ered as incompressible fiuid where cell motion is deter- 
mined by convective velocity field [1, 0, Hi • In the models 
of this type several convection equations describe dynam- 
ics of different tumour cell types or phases. Local changes 
in cell population lead to variations in the internal pres- 
sure, which in turn induces motion of tumour cells. These 
models also include reaction-diffusion equations for spa- 
tial nutrient and drug distribution within the tumour. 
However convection-domination models usually neglect 
tumour own cell motility. There are rather few models 
which take into account both convective and random tu- 
mour cell motion as well as chemo- or haptotaxis [ol, |l0| . 

For simulation of tumours with large random cell 
mot ility Fisher-like equations have usually been used 
prl . Il2| . In these models a logistic shape for the tumour 
cell proliferation rate is often used to prevent unlimited 
growth of local cell density. This rather artificial ap- 
proach does not take into account nutrient distribution 
inside the tumour and completly disregards the main can- 
cer property - unlimited growth in the case of a sufficient 
level of nutrients. 

With a rare exception the majority of papers on can- 
cer simulation consider solid tumours as compact objects 
with the total tumour cell density close to the maxi- 
mal possible cell density in the tissue. However there 
is another tumour type, namely infiltrative tumours, for 
instance glioma, characterized by a rather low value of 
tumour cell relative density and a large volume of pene- 
tration in normal tissue provided by high individual cell 
motility There are several mathematical models 

which describe tumours of this type (see, for example, 
[141]). However, spatio-temporal dynamics of nutrients is 
not taken into account there. 

In the present study a 1-D mathematical model for the 
infiltrative tumour is developed. The model exploits the 
fact that tumour cell division is possible only in the case 
of sufficient nutrient concentration. Thus the account of 
the nutrient spatial distribution effect on the tumour de- 
velopment is the central point of the model considered. 
The model consists of reaction-diffusion equations for the 
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tumour cell density and nutrient concentration. Depen- 
dence of the tumor front propagation rate on the model 
parameters is obtained both analytically and numerically. 
Properties of the propagating wave of the dividing tu- 
mour cells described by the model are compared to the 
solution of the Fisher equation. 



II. 



THE MODEL 



We consider a tumour as a colony of living and dead 
cells surrounded by normal tissue. Living cells divide 
with a constant rate and can move in a random way. 
In the case of nutrients lack tumour cells start to die. 
We take into account in the model that though a vari- 
ety of nutrients are necessary for the tumour growth, the 
oxygen shortage is mostly responsible for the cell death. 
We consider a tumour growing in a normal tissue with 
rather poor capillary system, so oxygen diffuse from the 
blood vessels located sufficiently away from the tumour. 
Oxygen consumption by normal cells which do not pro- 
liferate is neglected as dividing cells consume nutrients 
much faster than non-proliferating ones. Normal tissue 
surrounding the tumour is also supposed not to hinder 
neither cancer cell motion nor proliferation. We restrict 
our analysis to the case of a single spacial dimension. 
More formally the model is based on the following phys- 
ical and biological assumptions 

• A tumour grows as a colony of living and dead cells 
with densities a and m correspondingly. 

• Living cells divide with a constant rate B and in 
the case of low nutrient concentration s starve to 
death. 

• The death rate P{s) depends on the nutrient con- 
centration in a threshold manner which is described 
below. 

• Only random motility (diffusion) with a constant 
coefficient Da of tumour cells is considered in the 
model. 

• Oxygen is considered as a crucial nutrient for tu- 
mour growth. Its distribution is governed by dif- 
fusion Dg and consumption by the living tumour 
cells only according to linear law qa. 

Using these assumptions the governing equations can be 
written as follows 



da „ d^a , , „ 



dm 
'dt 



P{s)a, 



(1) 



ds d'^s 



where a, m and s are the proliferating cell density, the 
dead cell density and the nutrient concentration corre- 
spondingly. These equations are written already in a 
non-dimensional form. In order to estimate the corre- 
sponding parameters we take characteristic scales of time 
and length as To ~ 10^ s and L = 5-10^^ cm respectively. 
The cell density is scaled on amax = 10'' cells/cm"^ (max- 
imal cell density) and the concentration of oxygen in the 
tissue near blood vessels is supposed to be Smax = 10"'' 
mol/cm'^. In dimensional values the cell proliferation rate 
corresponds to the cell division frequency of the order of 
one division per 100 hours. The diffusion coefficients for 
oxygen and cells are equal to Dg — 2.5 • lO^^cm^/s and 
Da = 2.5 • 10^^ cm^/s respectively. Thus we obtain the 
following non-dimensional parameters of the problem 



Da - 10, 
Scrit =0.3, 



Ds = 10\ 
e = 0.01 



B = 0.1, 
g= 1.0, 



Prn = 0.2, 



(2) 



which will be referred to as a standard parameter set. 

The cell death rate is governed by P{s) which is a step- 
like function. For s > Scrit it is almost equal to zero and 
for s < Scrit it is greater than B. We will take it in the 
form 



Pis) = P„ 



1 - tanh[(s - Scrit)^-] 



(3) 



where Pm is the maximal value of P(s) and the parameter 
e defines the characteristic deviation of s form Scrit at 
which the death rate changes form the values close to 
zero to the values close to Pm- 

In Eq. ID) the second equation decouples from the 
rest set. The dead cells density profile determines the 
position of necrotic region inside the tumour. Therefore 
it does not affect the tumour spreading, which is governed 
only by the dynamics of proliferating cells and nutrient 
concentration and thus Eq. ([1]) is reduced to the following 
equations 



at 

ds 
dt 



Ba, 



(4) 



Ds 



' dx^ 



qa. 



At the avascular stage a tumour is supposed to have a 
spherically symmetric shape. However if the radius of the 
tumour is much greater than the characteristic scale, for 
which the distribution of cell density and nutrient con- 
centration change significantly, then a planar geometry 
can be considered. In view of our assumption that the 
tumour grows in the tissue where oxygen predominantly 
diffuses from the blood vessels located far away from the 
tumour, the set ^ can be solved in an infinite region. 
Thus we supplement Eq. ^ with the following boundary 
conditions on a and s 



0, 



for 



a = 0, 

S — Sn 



for 



(5) 
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Undoubtedly, any mathematical model of a biological 
system depends crucially on the choice of parameter val- 
ues. The parameters vary drastically depending on tu- 
mour type, localization, progression etc. The choice of 
parameters in the model is determined by our objective 
which is rather a qualitative description of an infiltrative 
tumour, but not a quantitative description of any specific 
tumor. Therefore the parameter values are taken in the 
experimentally observable range and they are not related 
to a certain kind of cancer. The other important limita- 
tion on parameters is that their values are chosen in such 
a way that the tumour cell density remains substantially 
smaller than the maximal value, which is typical for an 
infiltrative type of a tumour. 



III. TRAVELLING WAVE SOLUTIONS 

We seek the solution to Eq. ([4]) in the form of a wave 
travelling with a constant shape and speed c (i.e. an 
autowave). Then the governing equations ^ are reduced 
to a set of ordinary differential equations by introducing 
the automodel coordinate frame ^ — x — ct 

^ d^a da , . 



7^ d'^s ds 
with the following boundary conditions 



(6) 



0, 
a 



for ^ 



-oo. 



0, 
1 



for 



+00. 



(7) 

Here the constant a corresponds to the limiting constant 
value of the substrate concentration on — oo. We refer to 
the parameters c and a as internal parameters contrary 
to the control parameters, listed in On the first step 
we consider the asymptotic behavior of the solutions of 
Eq. ([6|) on ^ = ±oo. In order to do this we linearize Eq. 
([S]) near the values ([7]) which are stationary points of ^ 
and obtain two sets of ODEs with constant coefficients 



D. 



d'^a da 



+ c—-{B - P{cr))a- ^ 0, 



(8) 



^ d^s ds _ ^ 

U« c^TT qa = 0, 



and 



Da- 



,^Ba+ = 0, 



--+c—^qa =0. 



(9) 



According to linear differential calculus the solutions to 
these problems are sought in the form (a, a^, s, s^)'^ ~ 



exp(/i^^) which reduces the systems of linear differen- 
tial equations to eigenvalue problems for coefficients fj,^ 
as eigenvalues and constant vectors as eigenvectors. 
Eigenvalues fi^ indicate the rate of exponential conver- 
gence (divergence) of the solutions to the boundary val- 
ues ([7]) as ^ ^ ±oo and superscripts and '-' refer to 
the linearized problem on -|-oo and —oo respectively. 

For the case ^ — > — oo it appears that the linearized 
set of ODEs ([8]) has a singe solution with the rate of ex- 
ponential convergence to the boundary conditions given 
by A^r = (-C + ^yc-'+4.iPia)-B)Da)/2Da, two solu- 
tions unbounded for ^ — > — oo with the rates of exponen- 
tial divergence = (-c- y^c^ + 4(P(ct) - B)Da)/2Da, 
/ig = —c/Dg, and a single solution with a zero coefficient, 
/i4 = 0, of exponential growth. It can be shown that the 
stationary point {5*1 : a = 0, s = tr} of ([6]) is cither a 
saddle-node if is real positive, a stable node if fj,^ is 
real negative, or a stable focus if is a complex number. 

In the same manner, for the case ^ +oo we de- 
rive the set of ODEs ([9]) linearized near the bound- 
ary conditions ([7|) which has four linearly independent 
solutions bounded for ^ — s- -|-oo and characterized by 
the rates of exponential convergence /i^2 = (^c ± 

Vc2 - ABDa)/2Da, nt = -c/Ds, fit = 0. A simple 
consideration shows that the stationary point {S2 : a = 
0,3 = 1} is either a stable node or a stable focus. In 
the latter case a can become negative and therefore this 
solution is physically unrealistic. 

The eigenvectors associated with the lin- 
earized problems can be written as k^2 — 

(1, Mo.^M^s. Mo^M^z)'^: 



(0, 0, 



1, 0) . Here we have introduced a 



and kf = (0, 0, 

notation = \J -DsMi^2 + c/if 2- 

Taking into account that the solution of Eq. (|6I7|) exist 
only if the coefficient fi^ is a real positive number {SI is 
a saddle-node) and /x^2 3 ^"^^ real numbers {S2 is a stable 
node) we derive the following restrictions on the model 
parameters 



P(cr) > B, c> 2^/BDa- 



(10) 



The last condition implies that an automodel solution 
of Eq. ([4]) can propagate only with the velocity higher 
or equal to some minimal value Cmin = 2\fBDa- This 
property of the set ^ solutions is similar to the one 
of the KPP (Kolmogorov-Petrovskii-Piskunov) [l^ equa- 
tion and its special case. Fisher equation jlGj , which also 
exhibit autowaves. 

The first inequality in (IIO|) shows that since P(s) is 
a monotonic decreasing function there exist autowaves 
characterized by the residual concentration of the sub- 
strate left behind the travelling wave and this resid- 
ual concentration is less or equal to some critical value 
Omax = P^^{B), i.e. lim5__oos(f) = cr < Omax- The 
existence of both the parameter a and the limiting con- 
dition makes our model different from the KPP model. 

It can be also shown that for the travelling wave so- 
lution of Eq. ([5]) the substrate concentration profile is 
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a monotonic function of the coordinate ^. In order to 
do this we multiply the second equation in ([5]) by the 
integrating factor exp{c/ D^S,) and integrate it over ^ to 
obtain 



a(z)e"/^=^dz. (11) 



The integral in the right hand side of Eq. (fTT|) is always 
positive therefore is positive and s(^) is a monotonic 
growing function. 



IV. NUMERICAL SIMULATIONS 

We solve Eq. ([6]) subject to boundary conditions ([5|) 
numerically by using shooting and relaxation methods. 
Here we skip the description of these methods since they 
can be found in our previous papers (see [13 and refer- 
ences therein). The problem for numerical calculation is 
posed on a finite domain ^ G [— where L is taken 
to be sufficiently large, with a uniform grid on it. The 
integration step and L values are chosen in such way that 
both twice decreasing the integration step and twice in- 
creasing the integration interval results in a variation of 
the calculated parameters (such as wave speed) in the 
ninth significant digit only. 

The stability analysis of autowave solutions travelling 
with different velocities is studied numerically. Equations 
(l4|) are solved in a coordinate frame moving with a mini- 
mal speed (^ = a; — Cmint)- The boundary conditions ([5|) 
are imposed at the edge points of the space grid. For our 
numerical algorithm we use the method splitting with re- 
spect to physical processes. Initially we solve the set of 
ODEs which describes the cell birth and death processes 
as well as nutrient consumption by using the fourth or- 
der Runge-Kutta algorithm. As a next step, equations 
of mass transfer for oxygen and cells are solved with the 
Crank-Nicolson method of the second order approxima- 
tion in space and time. 

Typical solution profiles for the cell density, a(^), and 
a projection of phase trajectories into the plane (a, a^) 
are shown in the top and bottom figures [T] respectively 



for various values of the wave speed c > c„ 



where 



Cmin = 2-\/2 in this case. In the bottom figure [T] we also 
show the asymptotic behavior of the solutions of Eq. ^ 
[7]) for ^ ^ ±00 with straight lines which represent the 
trajectories of problems obtained form Eq. ([6]) by lin- 
earizing it near the boundary conditions ([7]). The cell 
density profile looks like a sharp peak so that a vanishes 
quickly form its maximum value in a relatively thin re- 
gion of the ^ coordinate. As we increase c the a(^) the 
profile becomes less localized and the maximum value of 
the cell density a(0) increases as well. It can be shown, 
however, that a(0) always remains less than one. The 
substrate concentration profile s(^) is a monotonic grow- 
ing function of coordinate ^ (not shown here for the sake 
of brevity), which is almost constant, s(^) ~ a for ^ < 
and exhibits a slow exponential growth to its boundary 



C = S 

c=10 
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FIG. 1: Cell density profile as a function of coordinate ^ in 
the moving frame (top figure) and projections of phase trajec- 
tories onto the plane (a, a^) (bottom figure) for standard pa- 
rameter values and various values of the wave speed c > Cmin ■ 



value (s = 1 as ^ ^ oo) for ^ > with the rate of growth 
c/Ds- The switching between these two regimes occurs 
in the thin region near C = where a reaches the maxi- 
mal value and coupling between the two equations in ^ 
is strong. 

In the top and bottom figures [2] we plot a typical so- 
lution profiles for the cell density a(^), and a projection 
of phase trajectories onto the plane (a, a^) respectively 
for the several values of speed below the minimal value 
c < Grain = 2 as shown in figure. In figure [5] it is clearly 
seen that at certain ^ values cell density becomes nega- 
tive which is not physically feasible. Therefore we con- 
sider the solutions for which a{£) > 0, travelling with the 
speed higher than or equal to a minimal value c,„i„. 

As a next step we investigate the dependence of the 
speed of the autowave solution as a function of param- 
eters of the problem. In the top figure [3] the depen- 
dence of c on (T is plotted for the physically feasible so- 
lutions travelling with the speed higher than Cmm (we 
refer to them as monotonic solutions) for standard pa- 
rameter values and several values of B as shown in Fig. 
[31 For fixed control parameter values there exists a fam- 
ily of solutions c((j) travelling with different speeds. For 
each value of B the speed is scaled on the minimal value 
Cmin = 2\/_B_Daj which is also shown with a dotted line 



marked "c 



All 



cease to exist when c/c„ 



monotonic solution branches 
I ratio drops down to 1. Be- 
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FIG. 2: Cell density profile as a function of coordinate ^ in 
the moving frame (top figure) and projections of phase tra- 
jectories onto a plane (a, a^) (bottom figure) for standard 
parameter values and various values of speed c < Cmin- 

low this value a family of the oscillatory solution (which 
are not physically feasible) branches emerges. These so- 
lutions are plotted in the bottom figure [3] with dashed 
lines. It is worthwhile noting that at c/cmin = 1 only the 
second condition in (fTO|l is violated, whereas the first con- 
dition is approached along the oscillatory branches when 
speed of the wave drops down to zero and a reaches max- 
imum which is shown in the bottom figure [3] 

We use the solutions obtained by solving Eq. ^ as 
the initial conditions for numerical integration of Eq. 

subject to boundary conditions (O in order to in- 
vestigate the stability of these autowaves. In Fig. [4] we 
show a spatial-temporal evolution of the initial profiles 
obtained by numerical integration of Eq. ([5]) with the 
speed equal to the minimal value, c = Cmm (top figure), 
and with the speed higher than minimum, c = 5 (bottom 
figure). The standard parameter values are used to per- 
form this calculation. We use the travelling coordinate 
frame, = x — Cmmt, so that the solution propagating 
with the minimal speed is presented in the top figure |4] as 
a standing wave, whereas the solution corresponding to 
c = 5 is shown in the bottom figure |3] as a wave travelling 
with the speed c = 5— c„im. The set ^ is integrated over 
the time periods of the order oi t — 4005""'^, where 
is the characteristic time scale for the instability growth 
in Eq. Results of numerical simulation of Eq. (jH) 

are shown in Fig. 3] for two values of the wave speed 



FIG. 3: Dependence of c/cmin on a for various values of B. 

and demonstrate that autowaves can propagate stably 
for long intervals of time without changing their speed 
and form. 



V. LIMIT e = 

In the limit e — > the governing equations ^ al- 
low the following simplifications: (i) the cell death rate 
function can be replaced by a Heaviside step function 
PmH{s — Scr), (ii) due to translational invariance of 
Eq. (|3]) we may require the travelling solution to sat- 
isfy the condition s(^)|^=o = Scrit- This implies that 
P(s(0) = for e < and P(s(0) = Pm for C > 0. In 
this case Eq. ([6]) can be considered in two semi-infinite 
intervals ^ G (— oo,0) and ^ e (0, oo) separately. This 
yields a set of ODEs similar to Eq. ([8]) for ^ < 0, where 
P{(j) is replaced by P,„, and a set (O for ^ > 0. The 
unknown functions a* and satisfy the boundary con- 
ditions ([7]) for ^ ±oo respectively and can be written 
explicitly as 

4 

u^(0=E"«^k±e'^?«, (12) 

i=l 

where u='= = {a^, a^", s^, s^)-^; and are taken 
from the section Hill (where P{<j) is replaced by Pm) and 
af are constant coefficients which are to be found. From 
boundary conditions ([7]) it follows that =0, =0, 
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FIG. 4: Dependence of cell density, a, on coordinate, ^, in 
the frame travelling with the minimal speed, and on time t 
for the standard parameter values. Initial profiles correspond 
to c = Cmin (top figure) c = 5 (bottom figure). 



a and a\ 



1. The other four constants can be 



found from continuity conditions across = for ^(0: 
a-(0) = a+(0), s(e): s-(0) = s+(0), and s^(0 = 5c„t: 
Sj (0) — st(0). This yields four equations 



0L-, = a, + a. 



2 ' 



^"i Ail = + + "3 mJ- 



(13) 



Our aim is to find the dependence of the speed of the 
autowave solution on parameters of the problem i.e. both 
the control parameters ^ and the internal parameter 
tr. The four unknown coefficients together with a 
make five unknowns. In order to find c we need one 
more equation which can be obtained by integration of 
the second equation from Eq. ^ over ^ with infinite 
limits and taking into account ([7]) and p^ . 



Mr 



"2 

/4 



(14) 



FIG. 5: Dependence of c/cmin on a for various values of B. 
The solid line represents the results obtained from Eq. psp 
and the dashed line shows the numerical data. 



The set of five equations (|13m4p is linear with respect to 



aj, ag , and a. It yields 



1 



(15) 



where fii is taken for ct = 0. The dependencies given 
by Eq. (fT5)) are presented in Fig. [5] with solid lines. The 
speed of the front is plotted vs a for various values of B. 
For each B the graph c{a) is scaled onto its critical value 
Cmin = 2\/DaB, which is shown in the figure with the 
dotted line marked as c/cmm = 1- The solutions travel- 
ling with the speed c > Cmin are described by Eq. (fT^ . 
For ^ > all coefficients are real and positive, there- 
fore a and s monotonically approach the limiting values 
([7]) as ^ tends to infinity. The solutions travelling with 
c < Cmin have complex coefficients /i^_3 with negative 
real parts, therefore for ^ > they exhibit oscillatory be- 
havior and for certain values of ^ > the cell density is 
negative contrary to reality . The dependence of c/cmin 
on (T is also shown in the same graph by the dashed line 
for the data obtained by numerically solving of Eq. ([B]) . 
The results obtained numerically and analytically are in 
good qualitative agreement. Quantitatively the influence 
of the finite length of the P(s) intermediate region (where 
it changes from P„i to zero) is stronger for small values 
of _B, as well as B close to P„i, whereas it is moderate 
for B close to Pm/2. 



VI. APPROACHING THE ASYMPTOTIC 
BEHAVIOUR 

It is known that for the Fisher equation the initially 
localized profile evolves with time in such a way that 
it approaches the automodel solution travelling with the 
minimal speed. Since the model considered in this paper 
has much in common with the Fisher equation, one may 
expect that it exhibit similar type of behaviour. Here we 
present the results of the investigation of the asymptotic 
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behaviour of the solution of Eq. (|4]) . The uniform distri- 
bution of the nutrient concentration, s = 1, and gaussian 
distribution of the hving cell density with the character- 
istic width A = 20 and amplitude 0.1 has been taken for 
the initial conditions (Fig. [6]). It has been shown that 
the process of the evolution of the initial profile can be 
split into two stages: initial and relaxation ones. 

The initial stage is shown in Fig. [5] (a), (c) and (e), 
where the maximal value of the cell density is plotted 
vs. t in figure (a) for the time interval [0, 100], the 
cell density and nutrient concentration profiles are shown 
for fixed moments of time: = (curves 1), t2 — 17 
(curves 2), ^3 = 30 (curves 3), ^4 = 55 (curves 4), and 
^5 = 100 (curves 5) in figures (c) and (e) respectively. 
Initially there is an excess of oxygen so that the increase 
of cell density due to division exceeds density decrease 
via both diffusion and cell death, which is negligible on 
this stage. Thus exponential growth of Umax is observed. 
Fast growth of the cell density is obviously accompanied 
with rapid consumption of nutrients in the regions where 
a is high. As a result s drops down below the critical 
value Scrit and the death rate P{s) hops from zero to its 
maximum value Pm in the regions where {Scrit — s) 3> e. 
This case is represented by the curve 3 in Fig. [6je) , where 
the critical value of s = Scrit is plotted by the dashed 
line. At this stage cells cannot reproduce themselves ef- 
fectively and cell density decreases due to high death rate 
which is becoming the dominating factor in the evolution 
of the a profile in the region close to amax- It changes 
the trend for amax from the exponential growth to practi- 
cally exponential decay with the index Pm — B resulting 
in the appearance of sharp peak in the amax(t) graph 
for t « 23. As a result the nutrient consumption de- 
creases, the nutrient concentration profile flattens due to 
diffusion from the regions, where it has not been con- 
sumed yet, and becomes monotonic. The cell density 
profile decays as is shown by curves 4 and 5 in Fig. EJe). 
For t €E [50, 100] there are oscillations in the amax value 
are observed due to redistribution of oxygen and cells to 
and from the region where a ~ amax- These oscillations 
vanish for t > 100, where the regime of relaxation to the 
automodel solution starts. The relaxation regime is char- 
acterized by a monotonic decay of amax to the maximal 
value of cell density for the automodel solution travelling 
with the minimal speed. This is illustrated in Fig. [5] (b), 
(d) and (f), where the dependence of amax on t (in figure 
(b)) and the solution profiles a{x,t) and s{x,t) (figures 
(d) and (f) correspondingly) are plotted for several mo- 
ments of time ti = 100, ^2 = 1000 and ^3 = 5000. In 
figure (d) we also plot by the dotted line the automodel 
solution profile travelling with the minimal speed. The 
dashed line in figure (f) corresponds to s = Scrit- Just 
as in case of the Fisher-KPP equation [l^ the relaxation 
of the cell density profile to automodel solution travel- 
ling with the minimal speed follows the power-mode law. 
However in our case the characteristic time of relaxation 
is proportional to Ds, what indicates that the dynamics 
of nutrient is important for adequate description of the 
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FIG. 6; Dependence of maximal cell density, Umax, on time, t, 
(figures (a) and (b)), the cell density profiles, a{x), (figures (c) 
and (d)) and the substrate concentration profiles, s{x), (fig- 
ures (e) and (f)) sampled at various moments of time. The cell 
density and nutrient concentration profiles are plotted in the 
coordinate frame travelling with Cmin- The dependencies a{x) 
are centered at the origin so that a(0) is equal to a maximal 
value of a over all x values for fixed t. 



tumour evolution. 

To complete the description of the tumour growth pro- 
cess, the full set of the governing equations (H]) was solved 
numerically. In Fig.[7]the overall cell density distribution 
a + m, which is comprised of both living a{x, t) and dead 
m{x,t) cells is shown. In this figure the active cell den- 
sity a is plotted by the dotted line and the overall cell 
density a -I- m - by the solid line. The profiles are sam- 
pled at fixed moments of time t = 50, 250, 500, 750, and 
1000. The initial conditions are similar to those, taken 
for calculations shown in Fig.[n]i.e. the substrate concen- 
tration is uniform in space and s{x, 0) — 1, the living cell 
density a is taken in a form of a gaussian distribution of 
the width A — 100 and amplitude 0.01 and the density 
of necrotic cells equals zero, m(x, 0) = 0. The dynam- 
ics of the gaussian profile a{x, t) and the corresponding 
dynamics of the substrate concentration are described 
earlier in this section. Therefore here we discuss how 
this dynamics correlates to the evolution of the overall 
tumour cell distribution and density of necrotic cells. At 
the initial stage while living cell density grows almost ex- 
ponentially with time the density of dead cells remains 
negligible. As tumour grows, the nutrient concentration 
drops below the critical value Scrit in the region where 
the majority of living cells are located. This results in 
the change of trend from exponential growth of living cell 
density to its almost exponential decay due to cell death. 
This process is obviously accompanied by fast accumula- 
tion of dead cells. The density m starts to grow rapidly 
so that a sharp peak in the m(x) distribution appears. In 
the course of this process most of living cells in the pri- 
mary tumour site die, forming the necrotic region shown 
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FIG. 7: The overall cell density concentration profiles a + m 
(solid line) and the living cell density profiles a (dotted line) 
vs coordinate x sampled at various moments of time: t = 50, 
250, 500, 725 and 1000. 

by a sharp peak in Fig. [T] 

A small fraction of living cells, that survived at this 
stage, spreads into surrounding tissue moving towards 
the source of nutrient. This results in the formation of a 
travelling wave of active cells on the tumour border. In 
the interior region a plateau of the dead cell m(x) density 
arises as a result of the living cell death due to the short- 
age of oxygen. As the profile of the living cells density 
approaches the automodel solution, the dead cell density 
behind the front tend to the stationary value rristat- In 
the case of e ^ this value is given by: 



mstat = ^ / (16) 

where a(^) is the living cells density profile, corre- 
sponding to the automodel solution of Eg. fl]) traveling 
with the minimal speed Cmin- 

The tumour cell density spatial distribution for t = 
500, 750 and 1000 is shown in Fig. [7l The total mahg- 
nant cell density a + m is close to the maximal possible 
value, equal to unity, near the tumour origin x — 500. 
It can be interpreted as the formation of a solid necrotic 
core in the primary tumour site. It is clearly seen that 
tumour grows only towards the nutrient source and out- 
side the primary site the overall malignant cell density is 
substantially less than unity. Thus tumour rather infil- 
trates neighbour normal tissue with a constant speed but 
not grows like a solid object. 

VII. DISCUSSION 

A mathematical model for the infiltrative tumour 
growth which includes distribution of living cells and oxy- 



gen in tissue is developed. The model adequately de- 
scribes a constant rate of the tumour linear size growth 
at the initial stage of neoplasm development and the for- 
mation of necrotic region in the tumour interior, observed 
in experiments [19]. 

In this model the existence of a family of autowave 
solutions travelling with different velocities is demon- 
strated. Their properties are investigated both numer- 
ically and analytically. It is shown that for fixed param- 
eter values the autowave solutions can propagate with 
speeds higher or equal to some critical minimal value 
Cmin- In this sense our model is similar to the KPP- 
Fisher model. However it has an essential new feature, 
compared to the latter, due to the presence of the second 
equation. Namely there is a close connection between 
the automodel wave velocity c and the residual nutrient 
concentration a behind the front. 

The evolution of the initially localized profile is inves- 
tigated. It is shown that initially localized distribution 
asymptotically approaches the automodel solution trav- 
eling with the minimal speed according to power-mode 
law. This again resembles the behavior of the autowave 
solutions of the Fisher equation. However in contrast 
to Fisher model in our case the rate of the convergence 
depends strongly on the parameters of the equation de- 
scribing the evolution of nutrient. From biophysical point 
of view the convergence of the initially localized profile 
to an autowave solution implies that if there is an initial 
center of tumor growth with nonzero density of active 
cells, it will eventually develop into a full scale neoplasm 
with a definite rate of linear growth of its size, accom- 
panied by formation of necrotic region inside of it. Our 
investigation shows that this process can be divided into 
three stages. The first stage is characterized by a fast 
growth of the tumour in the region with non-zero initial 
tumour cell density. This stage lingers until nutrient con- 
centration drops below the critical value. The first, fast 
growth stage is followed by second, intermediate, stage 
with a rapid death of living cells in the tumour original 
position and appearance of necrotic region there. The 
second stage lasts until the nutrient consumption by the 
tumour gets counterbalanced by its diffusion from the ex- 
ternal source. Finally, the travelling wave pattern which 
evolves to the automodel solution is formed and the ac- 
tive tumour cells start to spread with constant velocity 
towards the source of nutrients. Properties of the au- 
tomodel solution are basically determined by the equa- 
tion for the malignant living cell density and have much 
in common with automodel solution in the KPP-Fisher 
equation. It should be noted that in real biological sys- 
tems with finite size autowave regime can be unattain- 
able due to the lack of time and/or due to influence of 
the boundary. 

In the present study we focused our attention on tu- 
mours of infiltrative type, in which malignant cell density 
is substantially smaller than the maximum cell density 
in tissue. In this case it appeared possible to consider 
a relatively simple model, which takes into account only 
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individual cell motility and disregards convective fluxes 
arising due to cell division, essential for solid tumours. 
Actually, model simulations demonstrate that with the 
exception of the primary site of tumour growth the total 
malignant cell density is substantially smaller than the 
maximum cell density in tissue. What is more important, 
the density of dividing cells, which gives rise to convective 
fluxes in solid tumours, constitutes only several percent 
of the maximal tissue density. Obviously, with the other 
set of model parameters this condition may not be ful- 
fllled. Since the tumour cell distribution tends to the 
automodcl solution we were able first to investigate its 
properties and to determine whether the tumour belongs 



to the infiltrative type and thus the simplified model can 
be used, otherwise convective fluxes should be taken into 
consideration. 
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